o 



(N 



O 
q 

O 
• 1— I 

Oh. 



m 

00 
O 



Expectations from the Liouville-von Neumann Equation Using Chebyshev 
Expansion 

Giacomo Mazzi-'^'[3 and Benedict J. Leimkuhler-"^ 

School of Mathematics, University of Edinburgh, Edinburgh EH9 3JZ, UK 
(Dated: 26 March 2010) 

We consider a natural dimension reduction technique for the Liouvihe-von Neumann equation for a mixed 
quantum system based on evahiation of a trace formula combined with a direct expansion in modified Cheby- 
shev polynomials. This reduction is highly efficient and does not destroy any information. We demonstrate 
the practical application of the scheme with a model problem and compare with popular alternatives. This 
method can be applied to autonomous quantum problems where the desired outcome of quantum simulation 
is the expectation of an observable. 



I. THE LIOUVILLE VON-NEUMANN EQUATION 

An ensemble of N quantum systems each described 
by a wave function I^E*) may be expressed in terms of a 
density operator g defined by: 



N 



(1) 



The temporal evolution of this quantity is characterized 
by the Liouville-von Neumann equation: 



dg 
dt 



= --iiH, g], 



(2) 



where we have set h = 1. If is expanded using a (finite) 
approximate basis set . . . Iv^n}}, ([2]) may be viewed 

as an ordinary differential equation in a matrix argument. 
In the time independent case the solution for ^ may be 
written: 



-iHt ^ iHt 

goe 



(3) 



It is possible and sometimes preferable to rewrite ([3]) 
by introduction of the Liouvillian i = Id ® H — H (g) Id, 
where Id is the identity matrix, allowing us to recast g 
as a vector: 



(4) 



The main issue is then to evaluate the exponential of 
the matrix L. In the typical case, the size of g grows 
exponentially in the number of particles; at the same 
time, for a large system and using a typical basis set, H 
(and L) will be a structured sparse matrix. 

In the past decades many methods have been devel- 
oped for numerical evaluation of the matrix exponential^. 
However for large sparse matrices, the methods usu- 
ally applied are those based on expansion in Krylov 
subspace^Ti^. The main idea is to project onto the 
subspace: 

K„iiL,g) = spa.n{go, Lgo, L'^ go, . ■ . , L"^ go}- (5) 



To get a suitable basis set of ^ we may use the Lanc- 
zos algorithm^, as L is Hermitian. The Lanczos method 
is an iterative method, a very desirable feature in the 
context of large sparse matrices. For specific details see 
Ref. iQ and Ref. 0. However, because of the fact that 
all these methods involve the propagation of the matrix 
g, they suffer from requiring that matrix operations (or 
matrix-vector operations) be performed at each step of 
calculation. 

While the evolution of the system is described by the 
density matrix, the outputs we are interested in obtain- 
ing from quantum simulations are the expectations of ob- 
servables, these being the only quantities we can compare 
with the experiments. In the density matrix formalism 
the expectation value of an observable Q, associated with 
an operator Q is written as: 



(Q(t)) = Trace{^)(^)Q}. 



(6) 



^' Electronic Mail: IG.Mazzi@sms.ed.ac.ukl 



Whereas in general quantum simulations the equation of 
motion is solved for a quantity, g, which has dimension 
n X n, the types of outputs we are generally interested 
in are just one dimensional objects (jH]). In this paper we 
exploit this fact and design an algorithm that computes 
(almost) directly the evolution of the expectation value 
dl]), instead of the evolution of the density matrix ([3]). 
This approach does not lose any information of the orig- 
inal system (the only errors arise due to truncation), but 
at the same time the method provides a powerful com- 
putational tool, with potential dramatic reduction in the 
computational cost, especially when dealing with large 
matrices. The main idea of this approach is to exploit 
features of a Chebyshev expansion for the matrix expo- 
nential in (21 . It is important to remark that Chebyshev 
polynomials have been widely used to get polynomial ap- 
proximations of function of matrices, and especially for 
the matrix exponential in The terms of the Cheby- 
shev expansion can be constructed iteratively with the 
price being a sequence of matrix-matrix multiplications. 

Several methods to evaluate ^ are discussed in a re- 
cent monograph, see Ref. however our proposed Direct 
evaluation of the Expectation values via Chebyshev poly- 
nomial (DEC) method is different because it does not 
solve the evolution for the density matrix. Instead it ex- 
ploits the trace evaluation in ^ and with the evaluation 



2 



of just one Chebyshev expansion it allows the solution of 
at any time. DEC can be extremely powerful when we 
are only interested in the expectation values; if instead 
it is necessary to evaluate the evolution of the density 
matrix itself Q, then a traditional approach might be 
the best choice. 



If more than one observable is required it is still pos- 
sible to use DEC. The only difference with the single 
expectation case is that we need to store different sets of 
Tfc, one for each operator Q. 



A. Stopping Criterion 



II. THE DIRECT EXPECTATION VALUES VIA 
CHEBYSHEV 

The preliminary step of our method is to rescale the 
matrix within the interval [—1, 1], as outside this interval 
the Chebyshev polynomials grow rapidly, and the expan- 
sion becomes unstable; to do that we need to evaluate the 
two extremes of the spectrum of L. In order to obtain 
extreme values we propose, as already mentioned in the 
literatureifl', to perform a few steps of Lanczos iteration, 
as this provides a good approximation for the extreme 
eigenvalues, for small computational cost. If we define 
these two values as a and /3, i.e. /3 < 'j{L) < a, we may 
rewrite L as L = (5* Id — LgD), where D — {a — f3)/2, 
S = {a + /3)/2, and -1 < a{Ls) < 1. We may then ex- 
pand the exponential of Lg in the Chebyshev polynomials 
and we arrive at the following equation for g: 



iLt , 



£•0 



e-'*^ E Ck{tD)niLs)g„ , (7) 



with to = Dt. Both Ckito) and Tk{Ls) can be calculated 
itcratively: 

Ck{t) = {2-Skfi){-i)''Jk{t), (8a) 

Tk+i{x) = 2Tk{x)x - Tk-i{x), (8b) 

with initial values ro(a;) = Id, Ti(x) — x. Jk{t) is the 
fc-th Bessel function of the first kind. 
If we insert (O into ^ we find: 



{Q(t)) = Trace I 



k=0 



Ck{tD)Tk{L,)go \Q}. (9) 



By exploiting the linearity of the trace operation we can 
pull out of the trace all time dependent parts, and eval- 
uate once for all the coefficients Tk{Ls). In fact we may 
rewrite © as: 



m)) 



-its 



J2 Cfc(tD)Trace{rfc(L,)Q} 

fe=0 



(10) 



This is the key equation of the DEC method as it 
is possible to store an array of scalar values Tk — 
Trace{(Tfepo)Q}- All the time dependent terms are just 
scalar values that have to be multiplied by Tk to get the 
evolution of Q at any time: 



m)) 



-its 



k=0 



(11) 



The number of terms for the polynomial expansion 
in ([T]) depends on a prescribed tolerance e, and on 
the time tu- In other methods based on Chebyshev 
approximatioE^i^, the following has been suggested as a 
stopping criterion: 



S-t- ||c„„,,(tD)|| < S. 



(12) 



Due to the zeros of the Bessel function J, at fixed time tjj, 
(1121) may hold for some n, even though the expansion has 
not yet reached the convergence regime; it may happen 
that for ni > n we have that Cnj^ltn) > c„(i£)). To avoid 
this effect it is enough to use as a stopping criterion a 
combination of two Bessel functions; the cost of such a 
stopping criterion is that at most we need to perform an 
extra iteration step (j8ap . In our numerical tests we have 
used the following: 

Vl|Cn„„x-l(il5)P + l|Cnmax(*I?)^ll < £• 

(13) 

The total time r plays a role here, since the larger r the 
more terms (Tfc,Cfe) will be needed to get \ck\ below the 
threshold e. 



B. Computation of the Expansion 

In order to optimise the number of terms we evaluate, 
but without having to check at each step whether we 
have already evaluated enough terms Tk, we propose to 
evaluate first {Q{t)), at the final time r, and to store the 

^^max values of ffc. 

From equation (j8ap it is clear that Ck depends on the 
Bessel functions. If we look at the asymptotic behavior 
of the Bessel function of fist kind, for any A; G TV, we have 
that, for k fixedii: 



Jkit) 



1 



r(A: + 1) V2 



limt^O, (14) 



where r{t) is the Euler-F and for n € Z we have that 
r(n) = {n — 1)1. Equation (fT4|) shows that for any fc ^ 0, 
in a neighbourhood of t = 0, Jk{t) is increasing mono- 
tonically with respect to t. This behavior is maintained 
for the whole interval [0, j^] where is the first zero of 
the derivative of Jk ■ It is possible to show (see Ref. [Ill, 
Eq.9.5.2), that k < j'/.; consequently we can say that if 
(fT3| holds for a given Umax at r and t < rimax, then we 
are in the monotonically increasing region for Jnma.x and 
"Aimax+i- In this case, equation holds also for any 
t < r. 
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iteration Tk+i{Ls) = 2Tk{Ls) - Tu-i (I8b|. But what is 
actually needed in all our calculations is TkQo- Because 
of the linearity of the iterative expression, we may multi- 
ply To and Ti by go and then use (|8b[) directly on TkQo- 
The iterated operation in then just a matrix-vector mul- 
tiplication. 

After all the Tk needed have been stored, it is possible 
to get {Q{t)) at any time t e [0, r] by evaluating: 



{Q{t)) = 



-iDts 



ck{ts)fk 



(17) 



k=l 



where the Ck are evaluated iteratively via (|8ap , and rrimax 
satisfies (fT3|) for tZ). 

For the details of the algorithm we refer to the Appendix. 



FIG. 1. Example of few integer order Bessel Functions of the 
first kind. 

The Bessel Functions of the First Kind of integer order 
may be evaluated directly by using a three-term recur- 
rence relation 

2n 

Jn+l{t) = ~'-fri — Jn-1- (15) 

It is well known that (1151) becomes numerically unsta- 
ble for n > t, see Reflla To improve the method, 
we may exploit the linear nature of the iterative algo- 
rithm. It is possible to use Miller's algorithm, and to 
solve an inverted form of (fT5|) . i.e. to solve for J„_i 
given Jn, Jn+i—- When using Miller's Algorithm it is 
suggested to expand the number of terms (providing a 
sort of buffer), i.e. to start the backward iteration pro- 
cess from TTistart = u + r, where n is the actual order 
of the function we are interested on and r is some small 
expansion. In this case we need to know already from 
an a priori error analysis how many iterations need to be 
performed to to get below the threshold e. 

It is possible to prove that for the rescaled Hermitian 
matrix i^, when applied to a vector of unit Euclidian 
norm we have^: 

||P„_i(iL,)^.o-e-"^=eo|| < 4(ei-(*/2")^-^)™form > t 

2m 

(16) 

where P„j (t) is the order m expansion in Chebyshev poly- 
nomials. This equation indicates that there is a super- 
linear decay of the error when m > t. The formula 
may be derived by examining the asymptotic behavior 
of the Bessel functions (fT4|) . We may then use the rela- 
tion 4(exp{l — (T/2m)^} 2^)™ < e to approximate m. 

C. EfFicient Implementation 

The cost of DEC is all in the first step. Note that the 
cost of the evaluation of any itself is roughly equiva- 
lent to that of a matrix-matrix multiplication, as per the 



III. NUMERICAL EXPERIMENTS 

Nuclear Magnetic Resonance (NMR) is a spectroscopy 
technique that exploits the interaction between nuclear 
spins and electromagnetic fields in order to analyse the 
samples. The temporal evolution of such a system is 
described via a density matrix that has size {21 + 1/2)" 
where I is the spin and n the number of nuclei. The expo- 
nential growth of the size of g with respect to n impedes 
the use of simulations when dealing with systems involv- 
ing more than few (5-6) spins. Many attempts have been 
made to solve this (see e.g. Ref.'T^, RefflTl for recent ap- 
proaches), even using Chebyshev polynomials^^. These 
algorithms have been developed to simulate both liquid 
systems, where the Hamiltonian is generally time inde- 
pendent, and for Solid-State NMR. In the last case the 
Hamiltonian is time dependent due to the non averaging 
out of anisotropic interactions during the motion of the 
sample. When the Hamiltonian is time dependent it is 
not possible to apply DEC, as it is not possible anymore 
to isolate the time dependent part out of the trace. 

As in many other physical systems, nuclear spin dy- 
namics provides a perfect example to test DEC, because 
the final outcome of the simulations is an observable, the 
free induction decay (FID) signal, and this result is the 
sole important quantity, as it is the only data available 
from experiment. 

We have used DEC to evaluate this quantity: 

/(t) -Trace {e(i)/p} (18) 

g{t) can be written as combination of Pauli matrices, and 
Ip is the shift up operator: Ip = Ix + ily As Hamiltonian 
we assumed a sum of isotropic chemical shift and the 
isotropic term of a pair interaction called Homonuclear 
J-couplings, that depends on the inner product /-,- • Ij^: 

n n 

H = -Y.^,I]+Y.Jii^yIi (19) 

For the initial density matrix we set qq = —ly, that is the 
result of the application of a so called x-pulse to a sample 
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nz = 6112 nz = 180096 



FIG. 2. Left: Sparse structure of the Liouvillian (n = 1024, 
nz — 6112) for a system of 5 spins. Right; Structure of the 
Liouvillian (n = 16384, nz = 180096) for a system of 7 spins. 
Approximately each spin is interacting with half the other 
spins in both cases. 



already under the effect of a strong constant magnetic 
field along the z directior^i^. This is the usual initial 
condition when the acquisition of the signal starts. 

An illustration of the structure of the Liouvillian ma- 
trix is presented in Figl2] The sparsity depends on the 
number of interactions among the spins. In most cases 
the J-coupling interaction matrix J is relatively sparse. 
In our numerical simulations each spin interacts with 
about half of the other spins. 

Due to the fact that our implementation involves only 
matrix-vector multiplication, techniques developed both 
for structured and unstructured sparse matrices may be 
exploited. 

For comparison of computational costs we tested this 
method with an increasing number of spin particles us- 
ing different methods to evaluate the exponential. In 
particular to examine the error we compared DEC with 
the expm function of MATLAB, that uses a scaling and 
squaring algorithm with Fade' approximation. In this 
way we evaluate once for all U — e^*^''* where dt is the 
step size of the simulation, and then at each time-step 
we propagate g: 

Qn+l =Ugn (20) 

It is well known that in terms of computational costs 
this simplistic approach performs poorly, so we consid- 
ered two more realistic model reduction algorithms, both 
based on a Krylov subspace expansion. 



A. First Alternate Method: Lanczos Iteration 

The first one is well known in the literature^ii^. It evalu- 
ates, through a Lanczos algorithm, an orthonormal basis 
Vm of the Krylov subspace Km{L, Qo)- 

An approximation for g(t) is: 

g{t) = e-'^^'go « \\go\\V^e-'^"^'ei, (21) 



where Tm and Vm come from the Lanczos algorithm. 
The Lanczos algorithm provides an orthonormal basis set 
Vm for the Krylov subspace KmiL, go) via a three-term 
recursion^: 

Pj+iQj+i = Lqj - a-jq-i - PjQj-i qi = go- (22) 

Tm is a tridiagonal matrix of size mxm, and ei is the first 
vector of the canonical basis of size n. This technique is 
very powerful for short time simulations, because with 
few iterations m it is possible to have remarkably good 
approximations, but for longer times larger Krylov sub- 
spaces would be needed to stay close to the real solution. 
On the other hand if we do not consider enough terms in 
the Lanczos algorithm for longer times, ((2T|) is no longer 
a reliable approximation. 

It is possible to set a stopping criterion for the Lanczos 
iterations^; for a given t we can find m such that: 

t[Tm]m+l,7n\\e ||m,l < £• (23) 

To assure a good approximation throughout the whole 
simulation we set m to verify ([23l) ioT t = t the total 
time of the simulation. In this way we have the certainty 
that the Krylov approximation error is below e for all 
time t < T. The equation (|23|) involves the evaluation 
of the exponential of a tridiagonal matrix at each step 
of the Lanczos method; when m is large this operation 
may become a serious bottleneck for the whole simula- 
tion. However in our numerical tests we found out that 
the order of magnitude of the required m was roughly 
m ~ n/10, where n is the size of the Liouvillian. For 
this reason although ([23]) is the proper way to stop the 
iterations, due to the fact that we were more interested 
in costs comparison than in error analysis, we chose ar- 
bitrarily m = n/10. 

B. Second Alternate Method: Zero Track Elimination 

The second method considered for comparison is a 
technique recently developed especially for NMR simu- 
lation called "Zero Track Elimination " , (ZTE)12. This 
technique is based on the idea of pruning out the ele- 
ments of g{t) which do not belong to K{L, go). In order 
to reduce the steps needed to evolve the full system, we 
monitor the elements of g{t) that stay below a chosen 
threshold ^ during this first evolution steps and intro- 
duce structural zeros based on these observations. The 
evolution is then performed in this reduced state space 
{gz,Lz)- 

The idea is extremely appealing, as once the propa- 
gator for Lz is evaluated all the subsequent steps have 
the cost of a reduced matrix- vector, and it is possible to 
use (PO)) for the reduced system. It is claimedii that the 
error of such an approximation is similar to what would 
be obtained by considering in the Krylov expansion the 
contributions coming from high values of n in L^gij^. 
There are however some drawbacks: 
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matrix size 



FIG. 3. Logarithmic comparison of computational costs for 
Lanczos, Zero Track Pruning (ZTE) and Direct Expectiations 
via Chebyshev (DEC), with dt = 0.1 N = 1000. 



• for this method there is no available convergence 
theory; 



• the performance depends strongly from the initial 
condition qq, and on H. As expected, in our tests 
the size of the reduced system could change by a 
factor 2 depending on the number of interacting 
spins. 



C. Summary of Results 

The error-to-cost (measured in CPU time) diagrams 
are shown in Figure [3] for all the methods described. It 
is clear in this example that DEC is more than an order 
of magnitude more efficient than the alternatives. Obvi- 
ously especially when the system is over-reduced to slash 
the computational costs, paying a price in terms of error, 
DEC still mantains the error below the expected toler- 
ance. To avoid instabilities coming form the evaluation 
of the Bessel functions in these numerical tests we set the 
tolerance to be e = le — 7. 

DEC performs at its best for short time simulations 
(i.e. when total time r is small), so that we do not need 
to evaluate a large number of Tk , and when at the same 
time the use of small time step dt is required, as the cost 
for any step after the first is negligible. For instance, 
while for all the other methods the cost of a 1000 step 
simulation with = 0.1, is comparable with a simulation 
of 1000 steps with dt = 0.1, for DEC there is a gain of 
an order of magnitude, see Fig|31 
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-•^ Lanczos dt=0.01 




DEC dt=0.1 
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FIG. 4. Logarithmic comparison of computational costs for 
Lanczos and DEC when simulating for the same number of 
total steps TV but with different stepzise dt. N = 1000 in both 
the cases. 



TABLE I. Comparison of computational costs, 
dt = 0.1, N = 1000 



full size 


cxpni 


Lanczos 


Reduced_^ 


ZTE 


Reduced_^ 


DEC 


16 


0.12 


0.15 


4 


0.14 


8 


1.99 


64 


0.14 


0.16 


8 


0.16 


24 


2.00 


256 


0.31 


0.20 


25 


0.20 


64 


3.17 


1024 


4.27 


0.93 


102 


0.56 


160 


4.91 


4096 


191.4 


12.66 


409 


4.43 


432 


7.15 


16384 




309.02 


1638 


298.57 


3296 


22.34 


65536 












138.14 


262144 












1064.7 



Uize of g for the reduced system 



IV. CONCLUSION 

In this article we have presented a new method for sim- 
ulation of observable in a mixed quantum system. By 
expanding the exponential of the Hamiltonian in Cheby- 
shev polynomials, and exploiting the trace operation per- 
formed when evaluating the expectation value of an ob- 
servable, it is possible to reduce the evolution of any ob- 
servable to a one-dimension function that can be evalu- 
ated directly. 

We also presented an optimal algorithm to perform 
such a calculation, and show how this new method can 
easily compete in term of computational costs with a va- 
riety of model reduction approaches, whilst maintaining 
the approximation errors below a chosen threshold. 

CM. is very grateful to Arieh Iserles for useful sugges- 
tions at the starting of this work. 
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Appendix: The DEC algorithm 

We provide here a detailed description of the algo- 
rithm. 

inputs: L hermitian matrix n x n, go vector 
outputs: expectation value f{t) evaluated at f{jAt), j = 
1,...,N. 

1. evaluate a, [3 via Lanczos s.t. a = min(Aj), /3 = 

max(Ai); 

2. scale L and get Lg 

3. evaluate Tq = Id^O; = Lggo 

4. while ||cfe|| < s 

5. Ck = {2 — Skfl){—i)'' Jk{TS), T = total time 
6- Tk+i = 2TkLs - Tk-i 

7. store fk = Trace{{TkQo)Q} 

8. end 

9. for j = 1 : N 

10. re-evaluate the Ck at dilFerent time t = jdt 



11- f{3) = T.Ckfk 

12. end 
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